rm(list = ls())
library(ggplot2)
library("reshape2")
library("plyr")
library(shape)
#library(xlsx)
library(foreign)
library(grid)
library(ggpubr)

data <- read.dta("consolidated_elec_coal_final.dta")

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)

#1. Not-Interpolated Data

df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)

###########################################
# Scatter Plot - Total and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_total_noipo, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_noipo, ymax=ub_total_noipo),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_rural_noipo, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_noipo, ymax=ub_rural_noipo),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=elecrate_urban, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_urban_noipo, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_noipo, ymax=ub_urban_noipo),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c_noipo, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c_noipo, ymax=ub_total_c_noipo),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c_noipo, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c_noipo, ymax=ub_rural_c_noipo),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c_noipo, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c_noipo, ymax=ub_urban_c_noipo),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE Year FE) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c_noipo_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c_noipo_rcy, ymax=ub_total_c_noipo_rcy),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################


# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c_noipo_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c_noipo_rcy, ymax=ub_rural_c_noipo_rcy),alpha = 0.25)+ 
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_elecrate_urban!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c_noipo_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c_noipo_rcy, ymax=ub_urban_c_noipo_rcy),alpha = 0.25)+ 
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


#########################################
#########################################
#####LOWESS##############################
#########################################
#########################################

df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)

df <- subset(df, log_total_cumul!="NA")
df <- subset(df, elecrate_total!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$elecrate_total, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)

df <- subset(df, log_total_cumul!="NA")
df <- subset(df, elecrate_rural!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$elecrate_rural, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# Scatter Plot - urban and log cumul capa #
###########################################


# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)

df <- subset(df, log_total_cumul!="NA")
df <- subset(df, elecrate_urban!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$elecrate_urban, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=elecrate_urban, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_elecrate_total, nreps=100, confidence = 0.95,span=0.9)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_elecrate_rural, nreps=100, confidence = 0.95,span=0.9)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_elecrate_urban, nreps=100, confidence = 0.95,span=0.9)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE Year FE) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_elecrate_total, nreps=100, confidence = 0.95,span=0.1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################


# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_elecrate_rural, nreps=100, confidence = 0.95,span=0.1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_elecrate_urban!="NA")
df <- subset(df, rcy_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_elecrate_urban, nreps=100, confidence = 0.95,span=0.1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


#df <- subset(df, r_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

#2. Big Coal Reserve Countries

rm(list = ls())
library(ggplot2)
library("reshape2")
library("plyr")
library(shape)
#library(xlsx)
library(foreign)
library(grid)
library(ggpubr)

data <- read.dta("consolidated_elec_coal_final_bigcoal.dta")


# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)


#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

###########################################
# Scatter Plot - Total and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_total, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total, ymax=ub_total),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_rural, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural, ymax=ub_rural),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_urban, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_urban, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban, ymax=ub_urban),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
###########################################
###########################################
###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################
###########################################
###########################################
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c, ymax=ub_total_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c, ymax=ub_rural_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c, ymax=ub_urban_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE Year FE) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c_rcy, ymax=ub_total_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c_rcy, ymax=ub_rural_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_urban!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c_rcy, ymax=ub_urban_c_rcy),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# LOWESS Plot - Total and log cumul capa #
###########################################

library(spatialEco)

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_total!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_total, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_rural!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_rural, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_urban!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_urban, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_total, nreps=100, confidence = 0.95, span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "loess", formula=y~x, color="black", size=1,span=0.8)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country Year FE ) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_total, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p1 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p1


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_total_cumul!="NA")
df <- subset(df, rcy_ipo_elecrate_urban!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 45)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2


#3. Per Capita Capacity (logged)

rm(list = ls())
library(ggplot2)
library("reshape2")
library("plyr")
library(shape)
#library(xlsx)
library(foreign)
library(grid)
library(ggpubr)

data <- read.dta("consolidated_elec_coal_final_logcapa.dta")


# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)


#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

###########################################
# Scatter Plot - Total and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_total, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total, ymax=ub_total),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_rural, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural, ymax=ub_rural),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_urban, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_urban, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban, ymax=ub_urban),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
###########################################
###########################################
###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################
###########################################
###########################################
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c, ymax=ub_total_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c, ymax=ub_rural_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c, ymax=ub_urban_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE Year FE) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c_rcy, ymax=ub_total_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c_rcy, ymax=ub_rural_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_urban!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c_rcy, ymax=ub_urban_c_rcy),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# LOWESS Plot - Total and log cumul capa #
###########################################

library(spatialEco)

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_total!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_total, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_rural!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_rural, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_urban!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_urban, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_total, nreps=100, confidence = 0.95, span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "loess", formula=y~x, color="black", size=1,span=0.8)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country Year FE ) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_total, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p1 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p1


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_total_cumul!="NA")
df <- subset(df, rcy_ipo_elecrate_urban!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-0.5, 0.5)+
  ylim(-40, 45)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2

#4. By Region

rm(list = ls())
library(ggplot2)
library("reshape2")
library("plyr")
library(shape)
#library(xlsx)
library(foreign)
library(grid)
library(ggpubr)

data <- read.dta("consolidated_elec_coal_final_r1.dta")
#note: consolidated_elec_coal_final_r1 is the data for East Asia & Pacific
#use consolidated_elec_coal_final_r2 for Europe & Central Asia
#use  consolidated_elec_coal_final_r3 for Latin America & Caribbean
#use  consolidated_elec_coal_final_r4 for Middle East & North Africa
#use  consolidated_elec_coal_final_r5 for South Asia
#use  consolidated_elec_coal_final_r6 for Sub-Saharan Africa


# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)


#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

###########################################
# Scatter Plot - Total and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_total, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total, ymax=ub_total),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_rural, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural, ymax=ub_rural),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_urban, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_urban, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban, ymax=ub_urban),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
###########################################
###########################################
###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################
###########################################
###########################################
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c, ymax=ub_total_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c, ymax=ub_rural_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c, ymax=ub_urban_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE Year FE) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c_rcy, ymax=ub_total_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c_rcy, ymax=ub_rural_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_urban!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c_rcy, ymax=ub_urban_c_rcy),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# LOWESS Plot - Total and log cumul capa #
###########################################

library(spatialEco)

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_total!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_total, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_rural!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_rural, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_urban!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_urban, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_total, nreps=100, confidence = 0.95, span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "loess", formula=y~x, color="black", size=1,span=0.8)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country Year FE ) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_total, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p1 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p1


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_total_cumul!="NA")
df <- subset(df, rcy_ipo_elecrate_urban!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2




#5. By different time periods

rm(list = ls())
library(ggplot2)
library("reshape2")
library("plyr")
library(shape)
#library(xlsx)
library(foreign)
library(grid)
library(ggpubr)

data <- read.dta("consolidated_elec_coal_final_d1.dta")
#note: consolidated_elec_coal_final_d1 is the data for before 1980
#use consolidated_elec_coal_final_d2 for before 1980-2000
#use  consolidated_elec_coal_final_d3 for after 2000

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)


#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

###########################################
# Scatter Plot - Total and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_total, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total, ymax=ub_total),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_rural, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural, ymax=ub_rural),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
p <- ggplot(data = df, 
            aes(y=ipo_elecrate_urban, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(y=ihat_urban, x=log_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban, ymax=ub_urban),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
###########################################
###########################################
###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################
###########################################
###########################################
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c, ymax=ub_total_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c, ymax=ub_rural_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c, x=r_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c, ymax=ub_urban_c),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE Year FE) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_total_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_total_c_rcy, ymax=ub_total_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_rural_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_rural_c_rcy, ymax=ub_rural_c_rcy),alpha = 0.25)+   
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_urban!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

p <- ggplot(data = df, 
            aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(y=ihat_urban_c_rcy, x=rcy_total_cumul))+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(ymin=lb_urban_c_rcy, ymax=ub_urban_c_rcy),alpha = 0.25)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p

###########################################
# LOWESS Plot - Total and log cumul capa #
###########################################

library(spatialEco)

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_total!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_total, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_total, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Total)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - Rural and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_rural!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_rural, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Rural)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Scatter Plot - urban and log cumul capa #
###########################################
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, ipo_elecrate_urban!="NA")
df <- subset(df, log_total_cumul!="NA")

lo.df.f <- loess.boot(df$log_total_cumul, df$ipo_elecrate_urban, nreps=100, confidence = 0.95,span=1)
lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=ipo_elecrate_rural, x=log_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(0, 14)+
  ylim(0, 116)+
  scale_x_continuous(breaks=seq(0,14,2))+scale_y_continuous(breaks=c(0,20,40,60,80,100), limits=c(0,116))+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Logged)", y="Electrification (Urban)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_total!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_total, nreps=100, confidence = 0.95, span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_total, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "loess", formula=y~x, color="black", size=1,span=0.8)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p



###########################################
# Residual (Country FE only) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_rural!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_rural, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country FE only) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, r_ipo_elecrate_urban!="NA")
df <- subset(df, r_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$r_total_cumul, df$r_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.8)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p <- ggplot(data = df, 
            aes(y=r_ipo_elecrate_urban, x=r_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p


###########################################
# Residual (Country Year FE ) - Scatter Plot - total and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_total!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_total, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev


p1 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_total, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Total, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p1


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - rural and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_ipo_elecrate_rural!="NA")
df <- subset(df, rcy_total_cumul!="NA")

#df <- subset(df, r_ipo_elecrate_total!="NA" & r_total_cumul!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_rural, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_rural, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Rural, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2


###########################################
# Residual (Country FE + Year FE) - Scatter Plot - urban and log cumul capa #
###########################################

# Create Data frame
df<-data.frame(data)
df <- subset(df, region!="NA")
df$bmr <- as.factor(df$bmr)
df <- subset(df, rcy_total_cumul!="NA")
df <- subset(df, rcy_ipo_elecrate_urban!="NA")

lo.df.f <- loess.boot(df$rcy_total_cumul, df$rcy_ipo_elecrate_urban, nreps=100, confidence = 0.95,span=0.1)

lo.df<-data.frame(lo.df.f$fit)
lo.df$hi <- lo.df$y.fit +1.96*lo.df$stddev
lo.df$low <- lo.df$y.fit -1.96*lo.df$stddev

p2 <- ggplot(data = df, 
             aes(y=rcy_ipo_elecrate_urban, x=rcy_total_cumul))+
  geom_point(aes(color=bmr, shape=region), size=2)+
  xlim(-5, 5)+
  ylim(-40, 40)+
  geom_line(aes(x = x, y = y.fit), data=lo.df)+
  #geom_smooth(method = "lm", formula=y~x, color="black", size=1, se = FALSE)+  
  geom_ribbon(aes(x = x, y=y.fit, ymin=low, ymax=hi), alpha = 0.25, data=lo.df)+  
  labs(x="Cumulative Capacity (Residual)", y="Electrification (Urban, Residual)")+
  scale_color_manual(name="bmr", 
                     labels = c("Before 1980", 
                                "1980-2000", 
                                "Post 2000"), 
                     values = c("1"="#F8766D", 
                                "2"="#00BFC4", 
                                "3"="#7CAE00"))+
  guides(colour = guide_legend(title="Years", override.aes = list(size=4)))+
  guides(shape = guide_legend(title="Region", override.aes = list(size=4)))+
  theme_bw(base_size=10)+
  theme(
    axis.text.y = element_text(size=12, face="bold"),
    axis.text.x = element_text(size=12, face="bold"),
    #axis.ticks.x = element_text(size=12, face="bold"),
    #axis.line.x = element_line(color="black"),
    #axis.line.y = element_line(color="black"),
    axis.title=element_text(size=14),
    legend.text = element_text(size = 7, colour = "black"),
    legend.background = element_rect(fill = "white"),
    aspect.ratio = 0.7,
    #legend.position = c(0.13, 0.90),
    #legend.position="none",
    #legend.title = element_blank(),
    #panel.border = element_blank(),
    #panel.background = element_blank(),
    #panel.grid.major = element_blank(),
    #panel.grid.minor = element_blank()
  )

p2



